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It is possible in principle to probe the many-atom potential surface using density functional 
theory (DFT). This will allow us to apply DFT to the Hamiltonian formulation of atomic motion in 
monatomic liquids [Phys. Rev. E 56, 4179 (1997)]. For a monatomic system, analysis of the potential 
surface is facilitated by the random and symmetric classification of potential energy valleys. Since the 
random valleys are numerically dominant and uniform in their macroscopic potential properties, only 
a few quenches are necessary to establish these properties. Here we describe an efficient technique 
for doing this. Quenches are done from easily generated "stochastic" configurations, in which the 
nuclei are distributed uniformly within a constraint limiting the closeness of approach. For metallic 
Na with atomic pair potential interactions, it is shown that quenches from stochastic configurations 
and quenches from equilibrium liquid Molecular Dynamics (MD) configurations produce statistically 
identical distributions of the structural potential energy. Again for metallic Na, it is shown that DFT 
quenches from stochastic configurations provide the parameters which calibrate the Hamiltonian. A 
statistical mechanical analysis shows how the underlying potential properties can be extracted from 
the distributions found in quenches from stochastic configurations. 

PACS numbers: 05.70.Ce, 61.20.Gy, 71.15.Nc 



I. INTRODUCTION 

The potential energy surface underlying the motion of 
a monatomic liquid is composed of intersecting valleys 
in the many-atom configuration space. In Vibration- 
Transit (V-T) theory, the atomic motion is described by 
vibrations within a single valley, interspersed by tran- 
sits, which carry the system across intervalley intersec- 
tions [l], The pure vibrational motion is described 
by a tractable Hamiltonian which accounts for the domi- 
nant part of the liquid's thermodynamic properties. The 
remaining transit contribution is complicated but small. 
The vibrational Hamiltonian is calibrated by potential 
parameters evaluated at local minima (called structures) 
of the potential surface. These have been previously eval- 
uated for a model of liquid Na based on an interatomic 
pair potential [3|. Our goal now is to introduce first 
principles electronic structure calculations within density 
functional theory (DFT) to calculate the structural and 
vibrational parameters of the V-T Hamiltonian. 

In recent years, DFT has been used successfully in 
liquid studies to support the development of ab initio 
Molecular Dynamics (MD). The original Car-Parrinello 
method 0] was successfully applied to the melting of Si 
Q. Calculations of the MD trajectory on the adiabatic 
(electronic ground state) potential surface have beenper- 
formed for a broad spectrum of elemental liquids @, 0, § - 
This work has provided accurate andphysically revealing 
results for Al [9[, Fe [l(|, and Ge [TTJ] . The techniques 
we employ to calculate the supercell ground state energy 
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and Hellmann-Feynman forces are quite similar to those 
of Kresse et al. [12| ■ On the other hand, rather than 
proceed to MD calculations, our objective is to calibrate 
the (dominant) vibrational part of the liquid dynamics 
Hamiltonian from properties of local potential minima. 
We believe that the Hamiltonian formulation will use- 
fully complement ab initio MD in the study of dynamical 
properties of liquids. 

The procedure of quenching a system to its potential 
energy minima was introduced by Stillinger and Weber 
[HI, EL EH, and has become a valuable technique for 
studying systems with interatomic potentials. The tra- 
ditional method is to quench to many structures from an 
equilibrium MD trajectory, and then use statistical me- 
chanics to extract the underlying statistical properties of 
the potential surface [la LLZ|- I n V-T theory, however, 
we need very few structures — in principle only one for 
a given liquid at a given volume. Hence the traditional 
procedure for finding structures is inefficient, as a large 
number of MD iterations is required to bring the system 
initially to equilibrium as well as to avoid unwanted cor- 
relations between quenches. The problem is especially se- 
vere for DFT, where each iteration requires a converged 
total energy calculation, which is computationally very 
costly. 

We propose a simpler and more efficient method to 
probe the underlying potential landscape. Rather than 
quench from equilibrium MD configurations, we quench 
from configurations that are independent of interatomic 
interactions, and very fast to generate. Our purpose here 
is to demonstrate two properties of this quench method: 
that it is capable of producing the entire distribution of 
potential energy minima, and when used with DFT it can 
achieve our goal of ab initio calibration of the Hamilto- 



man. 

To perform the calibration, we will interpret structural 
data via the original hypothesis of V-T theory [l| : The 
potential energy valleys are classified as random and sym- 
metric. The random valleys numerically dominate the 
liquid statistical mechanics as N — > oc, and they all have 
the same potential energy properties as N — > oo; hence 
any one such valley may be used to calibrate the Hamil- 
tonian. The symmetric valleys have a broad range of po- 
tential energy properties, but make an insignificant con- 
tribution to the liquid statistical mechanics as N — ► oo. 
This hypothesis has been verified for the pair potential 
model of liquid Na at N = 500 0, and work in progress 
extends this verification to N = 4000 [H[ . 

The quench calculations are carried out for metallic Na 
at the density of the liquid at melt, using our model inter- 
atomic potential, and also using DFT. The Na potential 
was derived in pseudopotential perturbation theory, with 
an added Born-Mayer repulsion, and was calibrated from 
experimental crystal data at zero temperature and pres- 
sure The potential has since been shown to give 
excellent results for a broad range of experimental prop- 
erties of crystal and liquid Na (for a partial summary, 
see [13]). While there is no doubt DFT will provide ac- 
curate total energy results, we shall still need to verify 
that DFT quenches arrive at random structures, not sym- 
metric ones. This verification will be accomplished with 
the aid of the Na interatomic potential results. 

Our application of DFT to liquid dynamics theory is 
being pursued for a number of nearly-free-electron met- 
als and transition metals. A preliminary report has been 
presented on results for Na and Cu [2l|. We have not 
attempted to study elemental liquids whose equilibrium 
configurations are influenced by molecular bonding. Ex- 
amples are As, Se, and Te, which have strong and weak 
bonds, and Ge, whose liquid structure shows a contribu- 
tion from covalent crystal bonds. This anisotropic bond- 
ing will complicate the random structures underlying the 
motion in such liquids. This complication remains be- 
yond the scope of the present work. Moreover, since we 
have not yet presented structural data for other liquid 
metals, the present conclusions are strictly valid only for 
liquid Na. The results are expected to apply to many 
elemental liquids, perhaps all of them, but this extension 
is not demonstrated here. 

We consider a system of atoms in a cubic box with 
periodic boundary conditions. We construct configura- 
tions in which the nuclei are distributed uniformly over 
the box, within a constraint limiting the closeness of ap- 
proach of any pair. These are called constrained stochas- 
tic, or simply stochastic, configurations. The procedure 
is described in Sec. UH and the spatial uniformity is veri- 
fied by means of pair distribution functions. In Sec. IIII1 
our twofold purpose is addressed by two separate quench 
studies. In Sec. IIII At the Na pair potential is used to 
quench both equilibrium MD configurations and stochas- 
tic configurations. Comparison of the results will validate 
the use of stochastic configurations. In Sec. IIII Bl DFT is 
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FIG. 1: Mean and standard deviation of the pair distribution 
functions g(r) for 1000 stochastic configurations of N = 500 
atoms. The error bars indicate the standard deviation per 
histogram bin. The center of the nearest neighbor peak in 
g(r) for the quenched configuration of Na, shown in Fig. [5] is 
at n = 3.77 A. 



used to quench Na from stochastic configurations. Com- 
parison of the results with pair potential results will con- 
firm that the DFT structures are random and therefore 
calibrate the Hamiltonian. In Sec. HVl relations are de- 
rived between quenched distributions and the densities 
of states in the underlying potential energy surface. This 
analysis provides the statistical mechanical framework for 
interpreting results of the present quench technique. Our 
conclusions are summarized in Sec. IVl 



II. GENERATING STOCHASTIC 
CONFIGURATIONS 

Only minimal information is required to generate 
stochastic configurations: the number of atoms N and 
the system volume V, plus a distance of closest approach 
which is described below. Nothing of boundary condi- 
tions or the system potential has to be specified; how- 
ever, after the stochastic configuration is constructed, for 
all further calculations periodic boundary conditions are 
used. 

For a cubic cell with volume V , we construct a config- 
uration by choosing the particle coordinates at random 
over the cell. Randomness is important, as we shall use 
it in Sec. IIVI to determine the statistical weight factors 
for stochastic configurations. Next, a configuration is 
discarded if any two atoms are closer than a distance 
d. This is done for practical reasons: The self consis- 
tent field (SCF) calculation of DFT will not converge if 
atoms are too close to each other, and the pair potential 
at very small radii could lead to numerical instabilities 
in the conjugate gradient method due to the repulsive 
core. In practice, the excluded space can be very small. 
For Na we choose d = 0.4 A, so the relative excluded 
space (47r/3) d 3 /Va, where Va is the volume per atom, is 



FIG. 2: Mean and standard deviation of the pair distribution 
functions g(r) for 1000 stochastic configurations of N — 500 
atoms. The bins have a constant volume of Vb = (1/8) Va, 
where Va is the atomic volume of Na. The center of the 
nearest neighbor peak in g(r) for the quenched configuration, 
shown in Fig. [6] is at r\ — 3.77 A. For this particular binning 
volume, we have only 2 bins below the nearest neighbor peak. 



only 6.5 x 10~ 3 . Hence the stochastic configurations are 
expected to be spatially uniform to a very good approx- 
imation. 

To test the uniformity of stochastic configurations, 
we examine their pair distribution functions g(r). The 
conditional probability density g(r) is constructed as 
follows: Pick a system atom as central atom and de- 
note its position by r = 0. Make a set of bins la- 
beled b = 1, 2, • ■ • , in the form of concentric shells. Bin 
b has inner radius rj,, outer radius r&+i> and volume 
Vb = (47r/3) (rL_ x — rf). The pair distribution function 
g(r) has histogram rib (Va/Vj,), where rib is the number 
of atoms in bin b. Given the small size of d described 
in the previous paragraph, we expect g(r) to be nonzero 
even at relatively small radii. It is therefore important 
to normalize the bin count of bin b with the correct vol- 
ume of the bin, instead of using the approximate volume 
Airr^Arb as is often done. The bin contents are then av- 
eraged over the choice of each system atom as the central 
atom. While the bin radii are arbitrary, we usually take 
either Ar^ = rb+i — r& = constant or Vb — constant. 

For Na at N — 500, we constructed 1000 stochastic 
configurations and the g(r) histogram for each. With 
Ar& = constant, the mean and the standard deviation of 
the g(r) histogram are shown in Fig.[T] The blank space 
at small r is the empty sphere of radius d. The scat- 
ter at small r reflects the decreasing Vb as r decreases, 
and the corresponding decrease in rib- With Vb — con- 
stant, the mean and standard deviation of the g{r) his- 
togram is shown in Fig. [21 There the distance between 
points increases as r decreases, but the standard devia- 
tion remains nearly constant because rib remains nearly 
constant. The figures show the uniformity of stochastic 
configurations for r > d, with d being very small com- 
pared to the mean nearest neighbor distance. 



FIG. 3: Potential energy distribution of 1000 stochastic con- 
figurations of N = 500 Na atom systems. 

We also show the distribution of potential energies of 
the 1000 initial stochastic configurations. The mean of 
the distribution shown in Fig.[3]is at 3.75 eV/atom, which 
is considerably higher than the mean potential energy 
after the quench (see Eq. (Q])). 

III. PROPERTIES OF THE QUENCHED 
STRUCTURES 

A. Validation of quenches from stochastic 
configurations 

Figs. [Hand [5] compare two distributions of &o/N, each 
from 1000 pair potential quenches at N = 500. Fig. 2] is 
obtained by steepest-descent quenches from equilibrium 
MD at 800 K. This figure is an extension of the work 
reported by Fig. [5] is obtained by conjugate gradient 
quenches from stochastic configurations. We have veri- 
fied the equivalence of the two quench techniques for our 
system; for a related verification, see [22l [23l| . 

In each histogram we see the distinct random and sym- 
metric distributions, consistent with the V-T hypothesis. 
The random distribution is taken to be the dominant 
peak, out to where the histogram vanishes on either side. 
The symmetric structures are interpreted as the isolated 
parts of the histogram that are located outside the main 
peak on the low-energy side. 

The random distribution is numerically dominant and 
very narrow. The mean and standard deviation of each 
random distribution is given by 

, ( -183.29 ±0.50 meV/atom (MD) 
o/ ~~ \ -183.37 ±0.50 meV/atom (stochastic). 

Note that the dominant volume-dependent part of the 
pair potential is omitted here (see [24], Eq. (1.1) and 
Fig. 1). For this reason, most of the binding energy is 
missing from the pair potential energies in Eq. |T]). The 
mean value is the most accurate approximation to the 
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FIG. 4: Potential energy distribution of N — 500 Na atom 
systems after quenching from 1000 configurations drawn from 
a molecular dynamics trajectory at 800 K, or 2.16T m , using 
a steepest-descent quench method. 

thermodynamic limit value that we can make. The stan- 
dard deviation is the error expected from quenching only 
once and using that result. If N is increased, the mean 
value is expected to change slightly in converging to its 
thermodynamic limit, while the standard deviation goes 
to zero as N — -> oo. For a physical measure of the differ- 
ence in mean values, we note that the main contribution 
to the liquid thermal energy is the classical vibrational 
contribution 3fcsT, which is 95.90 meV/atom at T m . The 
difference in means is 0.08% of this. Experimental error 
in the thermal energy of elemental liquids at T m is typ- 
ically (0.1 - 0.5)%. Hence the two random distributions 
in Figs. [4] and [5] are identical to better than experimental 
error. 

Performing so many quenches has allowed us for the 
first time to see a clear and meaningful distribution of 
symmetric structures (compare for example @, H3, HH ) . 
Quenching from equilibrium MD yields 18 symmetric 
structures out of 1000 quenches, Fig. [U Quenching from 
stochastic configurations yields 23 symmetric structures 
in 1000 quenches, Fig. [5j Very approximately, the sym- 
metric distribution is constant and ranges from the bcc 
crystal, $ b cc /N = -196.12 meV/atom \2^, to the lower 
end of the random distribution. This broad distribution 
with few structures is consistent with the V-T hypothe- 
sis. If N is increased the symmetric distribution width is 
expected to remain the same, while the relative number 
of symmetries is expected to become negligible. In all 
these properties, the symmetric distributions in Figs. 0] 
and are the same to statistical accuracy. 

B. Calibration of the V-T Hamiltonian 

In order to calculate ab initio the thermodynamic 
properties of liquid Na for comparison with experiment, 
we quenched a stochastic configuration at Va — 41.27 A 3 
with DFT [2l| . The normal mode frequency spectrum 



FIG. 5: Potential energy distribution of N — 500 Na atom 
systems after quenching from 1000 stochastic configurations 
using a nonlinear conjugate gradient method. 

g(ui) was also calculated at this volume. For almost all 
monatomic liquids, the vibrational motion is nearly clas- 
sical at T > T rn . This means that the essential informa- 
tion required from g{uo) is the logarithmic moment of the 
frequencies, which provides the characteristic tempera- 
ture #o- Hence for calculation of thermodynamic prop- 
erties, the V-T Hamiltonian is calibrated from &o/N for 
the energy, and 9q for the entropy. Additional data which 
automatically accompanies the calculation of g{uj) will 
calibrate the Hamiltonian for nonequilibrium statistical 
mechanics. 

The DFT calculation is done with the VASP code 
using the projector augmented wave (PAW) method in 
the generalized gradient approximation (GGA) [2?], [28| . 
The planewave energy cutoff is 101.7 eV, the maximum 
core radius is 2.5 A, and the total energy convergence cri- 
terion is 10~ 8 eV. We use a r-centercd Monkhorst-Pack 
grid [2!| with 14 fc-points in the irreducible Brillouin 
zone. The total energy convergence for these parameters 
was carefully verified. Note that it is the large size of 
the real-space supercell which allows us to use few fc- 
points in comparison to the large number (several thou- 
sands) needed for crystal metal calculations with small 
unit cells. The quench is done by nonlinear conjugate 
gradient method. The system is considered quenched 
when the energy difference between subsequent config- 
urations is 10~ 7 eV or less. The DFT structure is at 
N = 150, a number large enough to get potential energy 
parameters accurate to a few percent, but small enough 
that convergence properties of the calculations can be 
studied. 

Because of the strong dominance of random valleys in 
the potential energy surface, the DFT structure is ex- 
pected to be random. To eliminate different zeros of en- 
ergy, we evaluate the energy difference 

A$ < c , (2) 

where the superscripts r and bcc represent respectively 
the random structure and the bcc crystal. The compari- 



__ J 12.75 meV/atom (pair potential) 



12.76 meV/atom (DFT). 



(3) 



The value for the pair potential is from the second mean 
in Eq. (JTJ), which is also calculated from stochastic con- 
figurations. The difference of 0.01 meV/atom is small 
compared to theoretical errors in A&o/N, and also com- 
pared to experimental error in the energy of liquid Na at 
melt. 

An independent confirmation is furnished by 9q. The 
Na pair potential value is from [3|]: 



98.4 K (pair potential) 
97.6 K (DFT). 



(4) 



The relative difference of 0.8% in 9q will make a corre- 
sponding difference of 0.3% in the theoretical entropy of 
liquid Na at melt. The difference is well within theoret- 
ical error in 9q, and is close to the experimental error in 
the entropy. 

The structural pair distribution G(r) is not a parame- 
ter of the V-T Hamiltonian, but G(r) has a role in density 
fluctuation phenomena, and it is therefore interesting to 
compare the DFT and pair potential results. The com- 
parison is shown in Fig. [SJ where the agreement is excel- 
lent. Notice the DFT curve (N = 150) has a small defi- 
ciency at the tip of the first peak, compared to the pair 
potential curve (N — 500). This deficiency is a small- 
N effect, and is observed also with the pair potential at 
N = 168, but not at N > 500 (0, Fig. 2). 



IV. EXTRACTING DENSITIES OF STATES 
FROM QUENCH RESULTS 

A. Quenches from equilibrium configurations 

In classical statistical mechanics, the partition func- 
tion for a single potential valley harmonically extended 
to infinity is e"' 3 * (T/9 ) 3N . The factor {T/6 f N ex- 
presses the vibrational motion. The transit contribution, 
which accounts for the valley-valley intersections, will be 
neglected here. The total liquid partition function Z is 



Z = 



G($ o ,0 o )e- /3 *° - Mod* 



T 



(5) 



where G(<E>0j #o) is the joint density of states for the col- 
lection of valleys. The normalization of G($>o,6o) is VV, 
the total number of valleys. The equilibrium statistical 
weight of a single valley is 



W eg {$ ,6 ) 



-/3*o 



(T/0o) 



:SJY 



(6) 



In equilibrium at T > T m , the probability of find- 
ing the system in d9 Q at 9q, and in d$o a t <&o, is 
P(§o,8 ) d0 o d$ Q , where 

P($o, e ) = G(* , 0o)W eq (<P o , 6 ). (7) 
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FIG. 6: (Color online) Structural pair distribution function 
G(r) for quenched structures from pair potential (solid curve, 
black, N = 500) and a single DFT calculation (broken curve, 
red, N = 150). 



Upon quenching from an equilibrium trajectory at T > 
T m , the structures sampled will exhibit a distribution 
proportional to P($o, Oo) & In view of Eqs. (JSJ) and © 
it follows that P(<&o,6*o) is insensitive to the normaliza- 
tion of G(&o,8 ). Therefore measurements of P(&q,9 ) 
can not be used to count the valleys. 

The probability of finding the system in d$o at $o is 
P(<E> ) d§ , where 



P($o) = J P($o,0o) d9 



(8) 



Upon quenching from an equilibrium trajectory at T > 
T m , the structures sampled will exhibit a distribution of 
$o proportional to P($o)- The distribution in Fig. [4] is 
proportional to P($o)- However, the density of states in 
$o is G($ ), given by 



G(*o) 



G($oA) d9 , 



(9) 



and differs from P( ( E > o) by the statistical weight 
W eq ($o, 9q). In our studies of liquid dynamics, the pur- 
pose of quenching is to determine the densities of states 
G($o) and G($o, 9q), because these are parameters of the 
Hamiltonian. These densities of states must be solved 
for from the observed distributions P($o) and P($o> #o) 
from the above equations. Even though the symmetric 
structures are supposed to be unimportant for the liquid 
as N — > oo, those structures will always be statistically 
important at T < T m , and will therefore be included in 
our analysis. 

Let us introduce subscripts r and s for random and 
symmetric, respectively, and write 



G(* o ,0o) = Gr(*o,0o) + G s ($ o ,0o), 



(10) 



and correspondingly Z = Z r + Z s . From Figs. 2] and 03 
P, (<I>o) at N = 500 has very small width, comparable 
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to experimental error in the internal energy of the liquid 
at melt, and this width is expected to decrease further 
as N increases [18j ]. These results suggest a model for 
G r (<I>o, Oo). Let us define the liquid $q as the mean $o 
for random structures when N — > oo, with a similar def- 
inition for 9 l . The model is 

GV($o, Oo) = Mr 8 ($ - $ ) 6 (9 - 6 l ) [l + O (N- a )] , 
G r ($ ) = M r 5 ($ - #{,) [l + O (i\T°)] , (11) 

where a > 0, and A/" r is the number of random valleys. 
From this it follows that Z r = M r e"^ (T/6 l ) 3N , and 
the random contributions to Eqs. © and ([5]) become 

5 ($ - $k) 5 (9 Q - 6k) r 



weights are 



1 + Z s /Z r 



^r($o) 



g (go - $b) 
1 + Z a /Z P 



[l + 0(A~ Q )] 



(12) 



Hence, to finite-A errors, the random valley Hamilto- 
nian parameters <£>q and 9 l are determined directly by 
quenches from equilibrium MD at T > T m . And, because 
of the form of Eq. (fTTj) for G r .($ , ^o), these observations 
will remain true when the statistical mechanics theory is 
improved to include transit effects. 

The symmetric density G s ( < i>o,6'o) apparently has TV- 
independent width with $o ranging from $g cc to <J>q. 
Symmetric structures with $o > $o exist, but they 
are rare for monatomic systems. The 9q dependence 
of G s ($o,^o) is not trivial. One expects that addi- 
tional (symmetry) parameters are important for symmet- 
ric structures. Nevertheless, G s ($o) and G s ($o,^o) are 
well defined, and can be extracted from quench data with 
the aid of Eqs. and ©. 



B. Quenches from stochastic configurations 



W r = 
W s = 



Vr 



M r V r +M s V s : 



MrVr+MsV s ' 

The probability distributions are 

P r {$0,9 ) = G r ($ ,9 )W r , 

P s ($o,0o) = G S (%,9 )W S . 



(13) 
(14) 



(15) 
(16) 



Hence the random and symmetric densities of states are 
each proportional to the distribution found in quenches 
from stochastic configurations. Applying the model of 
Eq. (EE) for G r ($ ,6>o) yields 

5 ($o - &„) 5 (6 - 9k) r 

(17) 

The conclusion here is the same as with equilibrium con- 
figurations, Eq. (fl"2"|) , that the parameters $g and 9 l are 
determined directly by quenches from stochastic config- 
urations, up to finite-A errors. The reason, of course, is 
the form of the model for G r ($ > 0o)> Eq. (fTTj) . For sym- 
metric structures, the above equations reveal two signif- 
icant points: 

1. Quenches from stochastic configurations can de- 
termine the magnitude of G s (&o,9q) relative to 
G r ($o, 9q), but only when W s /W r is known. 

2. The relation between G r ($0)#o) as determined by 
the two quench methods is unknown until W s is 
evaluated. 

These points are relevant to the distributions shown in 
Figs, m and 



On an equilibrium trajectory at T > T m , the probabil- 
ity the system is found in a given potential energy val- 
ley is W eq ($o, 9q) for that valley. The statistical weight 
is quite different for stochastic configurations. These 
configurations are uniformly distributed over configu- 
ration space, except for the small excluded Cartesian- 
space volume at each nucleus. Neglecting this constraint, 
the probability the system is found in a given potential 
valley is the valley volume divided by the entire 3A- 
dimensional volume. Let us denote the corresponding 
statistical weight factors as W r ($o,^o) an d W s ($o,0o) 
for random and symmetric valleys respectively. 

Because of their uniformity, the random valleys all have 
the same volume in the thermodynamic limit. To arrive 
at a complete solution, it is necessary to include the sym- 
metric valleys. Let us assume that they also have a uni- 
form configuration space volume. The number of random 
(symmetric) valleys is denoted M r (M s ), and the single- 
valley volume is V r (V s ). The statistical single- valley 



V. CONCLUSIONS 

In this article, we have addressed two main research 
goals: (1) The use of stochastic configurations to probe 
the distribution of potential energy minima, and (2) the 
calibration of the V-T Hamiltonian with DFT. 



A. Stochastic Configurations 

Quenches from stochastic configurations produce sta- 
tistically indistinguishable distributions of the potential 
energy, <f>o/N, compared to quenches from equilibrium 
liquid MD trajectories. We have demonstrated that 
quenching from stochastic configurations can be used 
to find the entire distribution of Na potential energy 
minima, i.e. they reliably reproduce the correct distri- 
bution of random and symmetric structures. This is 
illustrated by a comparison of $o/N distributions for 
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quenches from equilibrium MD and from stochastic con- 
figurations (Figs. H] and Sec. IIII A]) . Quenches from 
stochastic configurations yield an accurate random distri- 
bution for Na in agreement with the random distribution 
from quenches from equilibrium MD, Eq. (JTJ). 

Compared to generating and selecting configurations 
from equilibrium MD trajectories, our stochastic configu- 
ration method is significantly faster and more economical 
(Sec. QTJ). Stochastic configurations do not require inter- 
atomic potentials or costly equilibration and long simu- 
lation times. Simply generate random Cartesian coordi- 
nates for each atom under a minimal excluded-volumc 
constraint to eliminate particle overlap. Hence, this pro- 
cedure requires very little computational effort, an econ- 
omy that accommodates ab initio quench methods even 
for large systems. 

B. Calibration of the V-T Hamiltonian 

Calibration of the V-T Hamiltonian is based on the 
presumed dominance and uniformity of random valleys 
as N — > oo. This view is given mathematical expression 
in the model for G r (§ ,9 ), Eq. [[IT]). It follows that 
the thermodynamic limit parameters $q and 8 l are de- 
termined directly from data for either MD quenches or 
stochastic quenches, up to finite-iV errors. 

We have demonstrated that the DFT structure in 
Sec. IIII Bl is random by comparing the mean potential 
energy <I>o/iV with the pair potential results in Eq. |3]), 
the phonon moment 6q in Eq. (f5]), and the pair distri- 
bution function G(r) in Fig. O We conclude that being 



random, the DFT structure can provide ab initio calibra- 
tion of the V-T Hamiltonian. 

As verified by their pair distribution function, stochas- 
tic configurations have nuclei distributed nearly uni- 
formly over the system volume (Figs. [JJandEl Sec. PJ - 
Therefore among stochastic configurations, the statisti- 
cal weight for a many-atom potential energy valley is 
(nearly) proportional to the valley volume. This is in 
contrast to quenches from equilibrium MD, which require 
extensive modeling to extract the Boltzmann factor from 
the weight [16]. Given the statistical weight, characteris- 
tics of the underlying potential surface can be extracted 
from data acquired by stochastic quenches fSec. IIV[) . 

We do not suggest that DFT quenches from stochas- 
tic configurations will invariably arrive at random struc- 
tures, just as quenches from an equilibrium MD tra- 
jectory may result in a symmetric structure. Indeed, 
some symmetric structures have appeared in our DFT 
quenches. Precisely what is required to eliminate sym- 
metric structures from any collection of quench data is 
an ongoing research question. 
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